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ABSTRACT 

Evidence that the Universe may be close to the critical density, required for its expan¬ 
sion eventually to be halted, comes principally from dynamical studies of large-scale 
structure. These studies either use the observed peculiar velocity field of galaxies di¬ 
rectly, or indirectly by quantifying its anisotropic effect on galaxy clustering in redshift 
surveys. A potential difficulty with both such approaches is that the density parameter 
U 0 is obtained only in the combination (3 = if linear perturbation theory is 

used. The determination of the density parameter flo is therefore compromised by the 
lack of a good measurement of the bias parameter b , which relates the clustering of 
sample galaxies to the clustering of mass. 

In this paper, we develop an idea of Fry (1994), using second-order perturbation 
theory to investigate how to measure the bias parameter on large scales. The use 
of higher-order statistics allows the degeneracy between b and Sl 0 to be lifted, and 
an unambiguous determination of Uo then becomes possible. We apply a likelihood 
approach to the bispectrum, the three-point function in Fourier space. This paper is 
the first step in turning the idea into a practical proposition for redshift surveys, and 
is principally concerned with noise properties of the bispectrum, which are non-trivial. 
The calculation of the required bispectrum covariances involves the six-point function, 
including many noise terms, for which we have developed a generating functional 
approach which will be of value in calculating high-order statistics in general. 
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1 INTRODUCTION 

The assumption that bright galaxies are biased tracers of 
the mass distribution has featured strongly in theories of 
galaxy and structure formation in recent years. The idea 
that galaxies are biased even on large scale has been stud¬ 
ied for example by Peacock & Heavens (1985) and Bardeen 
et al. (1986), where galaxies are hypothesized as forming at 
high peaks of the density field filtered on relatively small 
scales. The concept of bias was then further extended to 
the extreme case where galaxies are painted on arbitrarily 
on the mass distribution. However, some constraints can be 
imposed, when the galaxy density is an arbitrary local func¬ 
tion of the mass density. In this case, the bias (here defined 
in terms of the two-point correlation function) must be a 
monotonic function of the spatial separation (Coles 1993) 
if it derives from a Gaussian field. Numerical experiments 
show that this also holds for the power spectrum (Mann, 
Peacock & Heavens 1997). 

It has been suggested that the bias mechanism can be 
modulated by environment dependent effects, as for exam¬ 
ple in the Cooperative Galaxy Formation model (Bower et 


al. 1993). The net result of this modification of the standard 
scheme is that the relationship between the density fluctua¬ 
tion field 5(x) and the galaxy fluctuation field 5 g (x) becomes 
non-local and the bias turns out to be scale-dependent. 
Other bias possibilities include dynamical friction effects, al¬ 
lowing galaxies to settle in clusters (Couchman & Carlberg 
1992). Note that even if the physical mechanism respon¬ 
sible for the bias operates only on relatively small scales, 
(e.g. within cluster cores), the bias parameter may differ 
from unity on much larger scales (Mann, Peacock & Heavens 
1997), with long wavelength modes having to have enhanced 
amplitude to fit the higher peaks of density in clusters. If the 
bias arises physically in the epoch of galaxy formation, it will 
evolve in time, and will approach unity if galaxy numbers 
are preserved (Fry 1996, Nusser & Davis 1994). This condi¬ 
tion may be broken by mergers (e.g. Matarrese et al. 1997), 
but also by luminosity evolution, if the sample is defined by 
a flux criterion, so it is safest to make no assumptions about 
bias evolution. 

There are clearly many ideas, but we have really very 
little idea of how bias works and evolves in practice, and 
there is therefore strong motivation to find ways to measure 
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it empirically from the galaxy distribution. The main moti¬ 
vation, however, comes from the desire to split the degen¬ 
eracy between S2o and b which arises in linear perturbation 
theory from dynamical studies of large-scale structure. Both 
direct studies of the peculiar velocity field (e.g. Dekel et al. 
1993) and studies of the distortion which peculiar velocities 
give to redshift-space galaxy maps (Kaiser 1987, Hamilton 
1992, Fisher, Scharf & Lahav 1994, Heavens & Taylor 1995, 
Ballinger, Heavens & Taylor 1995) yield a value of 


but not flo or b separately. The ignorance of the bias pa¬ 
rameter compromises all conclusions about the real prize - 
Qo, from large-scale structure studies. 

With good data, probably available in the next few 
years from the new generation of galaxy surveys, we show 
that the bias parameter can be measured by studying higher- 
order characteristics of the density field. More specifically 
this is done by analyzing the bispectrum: the Fourier trans¬ 
form of the three-point correlation function (e.g. Frieman & 
Gaztanaga 1994). This statistic has the advantage that the 
determinations of the bispectrum can be made uncorrelated 
on different scales, or if not, their correlation properties can 
be calculated, and the choice of modes to analyse can be 
made to maximize signal-to-noise. These allow the bias pa¬ 
rameter to be estimated via likelihood methods, importantly 
allowing error bars to be placed on the parameter estimates. 
The higher-order statistics work by exploiting the fact that 
gravitational instability skews the density field to high den¬ 
sities as it evolves. This behaviour can also be mimicked by 
non-linear biasing, but these two effects can be separated by 
the use of shape information: essentially non-linear biasing 
of a truly Gaussian field will lead to different shaped struc¬ 
tures from a non-Gaussian field arising from gravitational 
instability. 

It is worth making a few remarks about the advantages 
of performing this sort of analysis in Fourier space rather 
than via high-order correlation functions in real space. The 
main advantage, which is shared by the power spectrum, 
is that the estimates of the correlation functions in Fourier 
space can be made uncorrelated (or at least their correla¬ 
tion properties can be readily calculated), with enormous 
benefits for estimation of parameters and error bars. The 
second advantage is that the real-space determinations of 
high-order correlations are really only known well in the 
non-linear regime, where perturbation theory may not ap¬ 
ply. In Fourier space, there is a clear separation of scales 
where perturbation theory works and breaks down. 

Translating this idea into a practical algorithm requires 
a number of issues to be tackled, including the effects of 
a varying selection function and redshift distortions, and 
the non-trivial issue of what shot noise terms appear in 
the required six-point function. This paper is principally 
concerned with the error analysis; a subsequent paper will 
present a study of the other issues. 

The plan of this paper is as follows. In Section 2 we 
introduce the second-order perturbation expansion for the 
(unknown) matter density field and the (observable) galaxy 
density field. The necessary mathematical techniques to cal¬ 
culate the IV-point functions, including shot noise, are given 
in Section 3. Sections 4 discusses optimal methods and prac¬ 


tical implementation to a numerical simulation of the galaxy 
distribution, and Section 5 discusses the results and the pos¬ 
sibilities of practical application to real data. 


2 SECOND ORDER APPROXIMATION 

As Fry (1994) pointed out, since the degeneracy between S2o 
and b is an intrinsic feature of the linear theory, one needs 
to go to second order in order to separate the parameters. 

Under the assumption that the initial fluctuation are 
Gaussian and that structure grows by gravitational insta¬ 
bility, the three-point correlation function is intrinsically 
a second-order quantity, and should be detectable in the 
mildly non-linear regime where second-order perturbation 
theory is a reliable description. 

The data we will use to constrain the bias parameter 
will be the real part of the three-point function in fc-space: 

D a = R e (^fci<5fc2^fc3)’ (1) 

for closed triangles, a is shorthand for the triangle specified 
by some triplet of k-vectors. Sj: = J d 3 f 8(f) exp(—ik ■ r) 
is the Fourier transform of the fractional overdensity field 
8(x ) = n(x) / h — 1, where n is the number density of galaxies, 
and n its mean. D a is related to the bispectrum B(ki, k 2 , k 3 ) 
by 

( 5 ki 5 kjk 3 ) = (2n) 3 B(k 1 ,k 2 ,k 3 )8 D (k 1 +k 2 + k 3 ) = (D a ),(2) 

where 8 D (k) is the three-dimensional Dirac delta function 
and the () indicates ensemble average, or, by the ergodic 
theorem, the average over a large volume. In what follows, 
we may occasionally refer loosely to D a as the bispectrum. 
We expand the density field to second order as: 

8(x) = 8^(x) + S^ 2 \x), (3) 

where 8 ^ is 0(8 R-* 2 ) and represents departures from Gaus¬ 
sian behaviour. 

Note that what we are dealing here is not a mildly non¬ 
linear field, but a highly nonlinear field which is filtered 
on some smoothing scale. The operations of smoothing and 
gravitational evolution do not commute, so there is some 
possibility that perturbation theory might be inaccurate. 
However, this problem is not unique to second-order per¬ 
turbation theory: it exists in even for the evolution of the 
power spectrum, which for practical cases can be treated 
well by perturbation theory. Experiments on the bispectrum 
by others and ourselves (see section 4) show that perturba¬ 
tion theory for the implicitly filtered fields works well if the 
smoothed field is not too nonlinear. 

2.1 Unbiased case 

For an unbiased distribution (or for the matter distribution) 
we can expand (di^ife) to second order, obtaining: 

(8i8 2 8 3 ) ~ (<5‘ 1) ^ 1) 4 1) > 

+( < 5i 1) ^2 1 ^3 2 ^} + cyclic terms( 231, 312). 

(4) 

Only the second term of this survives. Since we will work in 
Fourier space we quote directly the expression for 8^?\ the 
Fourier transform of S^(x) (Catelan et al. 1995): 
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Figure 1 . The weak dependence of J on Qo for equilateral tri¬ 
angles. 



where 


J(ki, k 2 , fio) 


d 3 k a d 3 k b S D {k a + k b - k)J(k a , k b )5^5^\ 


(5) 


— i — s(n 0 ) + 


ki ■ k 2 
2k\k 2 


(ki k2_ 
\k 2 + fei 


) +s( n 0 ) 


( kik 2 \ 

^ kik 2 J ' 

( 6 ) 


B is a function which is 2/7 for an Einstein-de Sitter uni¬ 
verse, in which case the expression above reduces to the form 
originally obtained by Goroff et al. (1986). The useful fea¬ 
ture of this is that J is almost independent of flo, as seen 
in Fig. 1, so we can use the Einstein-de Sitter value with 
minimal error. 

The bispectrum for the unbiased case, in the absence of 
shot noise, is therefore 

( S ki S k 2 S ks) = (2tt) 3 [2J(ki,k 2 )P(ki)P(k 2 ) +cyc.(23, 13)] x 

5 D (ki + k>2 + fc 3 ), (7) 


where the linear power spectrum is defined in terms of the 
orthogonality of the transform coefficients, arising from ho¬ 
mogeneity: 

(S- kl 5- k2 ) = (27r) 3 P(k)5 D (k 1 + k 2 ). (8) 

For practical cases, the transform is made in a finite box, 
when the Dirac delta function is modified (see Section 3.3). 


2.2 Biased case 

We assume that the biased (galaxy) density is a local func¬ 
tion of the unbiased (matter) density. Galaxy properties will 
be identified by the subscript <?; matter fields have no sub¬ 
script. We make a Taylor expansion of the galaxy density to 
second order in so for consistency we must include an 
extra quadratic bias term b 2 : 

5 a{x) = f[S(x)] 


~ b\5{x) + ^b 2 S 2 (x) 

~ bid^\x) -|- biS (2 \x) -|- ^b 2 S (1)2 (x). 

(9) 

The linear bias parameter b i is to be identified with b. Notice 
that the expression above is not correctly normalised: a term 
b 2 {S^ 2 )/2 needs to be subtracted so that the biased field 
has zero mean. As we will be working in Fourier space, this 
extra term is irrelevant for all except the k = 0 mode, and 
will be ignored!/]. 

Following the same procedure as in the unbiased case 
we obtain to second order 

(S g iS g 2 S g3 ) = bUs^s^s^) 

+ 6? + eye. 

+ ^(^ 1) 4 1) 4 1)2 ) + cyc. 

( 10 ) 

Once again, the first term vanishes. We can re-write the last 
term by exploiting the cumulant expansion theorem for the 
four-point function obtaining the non-vanishing terms 

(sWx-^ 2 ) + 2(4 1) 4 1) ><4 1) 4 1) >- (ii) 

Therefore to second order the galaxy bispectrum is: 
( S gk 1 S gk2 S gk 3 '> = 

(2tt) 3 {6? [P(fci)P(fe)2J(fci, fc 2 ) + eye] 

+ b\b 2 [P(fci)P(fc 2 ) + eye] } 5°(k\ + k 2 + k 3 ). 

( 12 ) 

From this we see how the bispectrum may be used to esti¬ 
mate the bias parameters. J is almost independent of flo, so 
a set of triangles in k-space can in principle allow us to mea¬ 
sure the bias parameters 6i and b 2 . It is important to realise 
that it is necessary to consider triangles of different shapes 
to lift a new degeneracy between b i and b 2 , since triangles 
of a common shape have the same J and © then gives in¬ 
formation only on the combination b\(2b\J + b 2 ). The same 
sort of degeneracy has been noted for the 3-point function 
in real space (see Frieman & Gaztanaga 1994). The strong 
dependence on shape has been illustrated (e.g. Jing, Borner 
& Valdarnini 1995), but there appears to be little work on 
using different shapes to lift the degeneracy. We return to 
the physical reasons for having to take more than one shape 
in the discussion. 

We wish to estimate 6i and b 2 by maximum likelihood, 
so we need an estimate of the shot noise, and the covariance 
between bispectrum estimates. We turn to this in the next 
section, but here note that the bispectrum (strictly defined 
in terms of the mean) is a real quantity. This is a result of 
isotropy - the bispectrum does not depend on which direc¬ 
tion round the triangle is taken. Accordingly we use only 

* It is not obvious that a physical mechanism will lead to a func¬ 
tional dependence which can be expanded in a Taylor series. Ex¬ 
amples are the ‘weighted bias’ model studied by Catelan et al. 
(1994) and the ‘censoring bias’ studied by Mann, Peacock &; Heav¬ 
ens (1997). It is an open question whether such schemes might be 
approximated by a continuous S g (5) relation when smoothed on 
larger scales. 
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the real parts of ^ g %^ g % 2 ^ g % 3 as our data, as the imaginary 
parts just introduce more noise. 

Since P{k) is unknown, we write (^|) in terms of po¬ 
tentially observable quantities, the linear galaxy power spec¬ 
trum P g = b\P (equal to the observed power spectrum at 
adequate order for second-order perturbation theory).: 


( 5 gki 6 gk 2 5 gkJ ~ 

(2n) 3 {ci [2P g {ki)P g {k 2 )J{ki, k 2 ) + eye.] 

+c 2 [P g (ki)P g (k 2 ) + eye.] } S D (ki + k 2 + k 3 ), 

where 

1 _ b 2 

ci = —; c 2 = -j. 
b 1 b{ 


(13) 


(14) 


Under the assumption of a uniform prior for the quantities 
ci & C 2 , the a posteriori probability for the parameters is 
the likelihood: 


£(ci,C 2 ) = 


(27r)"2-(detC')s 

1 \ ' 


exp 


2 

ot/3 


(D a y a )C a/3 (Dp Vp) 


(15) 


where the means of the M data are y a {bi,b 2 ) = {D a ) and 
Cctf. 3 = {{D a — /j, a ){Dp — up)) is the covariance matrix, which 
is also a function of the bias parameters. The use of a Gaus¬ 
sian likelihood function is justified if there are many con¬ 
tributing data, whether or not they themselves have Gaus¬ 
sian distributions, provided the sensitivity to the param¬ 
eters is through the mean rather than the variance. It is 
straightforward to show that, by using the mean of many 
observed data drawn from the same distribution, the like¬ 
lihood is identical (by the Central Limit Theorem) as the 
Gaussian likelihood for the individual data, with the same 
variance and mean. In our case the variance is only weakly 
dependent on the parameters, through the shot noise, so our 
use of a Gaussian likelihood function is quite justified. 

Note that our data are real D a = Re{Z a ) = 
and we can write the covariance matrix for 
D a in terms of the complex Z a . 

{D a Dp) = \Re [{ZaZp) + (Z a Zp)] . (16) 

It should be obvious that we require the six-point function 
to calculate the covariance matrix, and we also need to in¬ 
clude shot noise in the bispectrum estimates. The former 
is sufficiently complicated by the presence of a very large 
number of noise terms that we present a general method for 
doing calculations of this nature in the next section. 


3 STATISTICS OF RANDOM FIELDS AND 
DERIVATION OF THE PENTASPECTRUM 


for example, from a flux-limited catalogue). These distribu¬ 
tions are required to calculate the covariance properties of 
the bispectrum. A reader interested principally in the appli¬ 
cation of the bispectrum may omit this at a first sitting. 

The complete description of a random field is given by 
the probability distribution functional P, which should be 
thought as the continuous limit of the joint probability dis¬ 
tribution in N points, as N —> oo. To describe the statistics 
of a general random field f{x), we need to know the func¬ 
tional V[f{x)], which specifies the distribution of the pos¬ 
sible values of / at each point. By statistical homogeneity 
P[f{x)] will be independent of x. Consider then the par¬ 
tition functional or generating functional of the correlation 
functions Z defined as a functional integral, as follows: 


Z[J{x)} 


D[f\P[f\ exp (i / d 3 xj{x)f{x)] 



i 


d 3 xj{x)f{x) 


(17) 


where J{x) is an external source perturbing the under¬ 
lying statistics and D[f] is a suitable measure such that 
the probability distribution is correctly normalized to 1: 
JV[f{x)]V[f{x)] = l. 

From this function we can recover the correlation func¬ 
tions p.N of the distribution and the connected correlation 
functions or kat as follows: 


• From the generating functional of the correlation func¬ 
tions 

y N = {fix ,)... f(x N )) = i~ N (18) 

evaluated at J = 0. 


• From the generating functional of the connected (or ‘re¬ 
duced’) correlation functions 

JC[J] = lnZ[J], (19) 


Kn {X 1 . . . Xn) — {f (^l) • • • f {xn )) connected 
-AT 5 n K.[J] 

1 sji... SJ N ’ 

( 20 ) 

evaluated at J = 0. 


From the above definitions it follows that a complete char¬ 
acterization of the statistics, i.e. the complete knowledge of 
the partition functional, requires the knowledge of the cor¬ 
relation functions of all orders. In fact it is usually possible 
to recover the partition functional as 

°° i n f 

Z[J] = 1 + ^— / dxl... dx n p,n{x\,. . . ,X„) X 

n=1 

J{xi) ■ ■ ■ J{x n ) (21) 

and similarly for KL[J\ in terms of the connected correlation 
functions, 


This section presents a mathematical method for calculat¬ 
ing A-point distributions in Fourier space for continuous 
and discrete fields, and deal also with the realistic case of a 
galaxy survey with a non-uniform number density (arising, 


K.[J] = ^2 ~[ / dx, .. .dx„K n {xi,... ,x„) x 


J{x\) ■ ■ ■ J{x n ). 


( 22 ) 
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The expansion of correlation functions of order N in 
terms of connected correlation functions of order M<N, i.e. 
the Cumulant (or Cluster) Expansion (e.g. Kendall & Stew¬ 
art 1977, van Kampen 1992), follows trivially from equations 
(21) and (22). The cluster expansion is widely used in sta¬ 
tistical mechanics where it constitutes the basis for approx¬ 
imation schemes in studying slightly non-ideal many-body 
systems like gases or fluids (Rice & Gray 1965). A review of 
the use of functional integration in cosmology can be found 
e.g. in Fry (1984), Bertschinger (1992), Matsubara (1995). 
It is worthwhile, at this point, to make some remarks on ap¬ 
plicability. In dealing with slightly imperfect fluids it is as¬ 
sumed that the interaction falls off faster than , this allows 
one to neglect the connected higher-order terms; in the case 
of gravity we must be aware that this is not generally true. 
The connected correlation functions kn, which are closely 
related to the connected Green’s functions of quantum field 
theory (e.g. Ramond 1989), are defined by the above proce¬ 
dure in such a way that any separable contribution of the 
type (f ){f N ~) has been removed. As a consequence, if 
any subset of the N points is set to infinite separation the 
connected correlation function of order N goes to zero. These 
connected correlation functions of order N for the density 
fluctuation field are often simply called ‘correlation func¬ 
tions’ in the cosmological literature (e.g. Peebles 1980). An 
important property of these functions is that if the under¬ 
lying statistics is Gaussian all the connected correlations of 
order higher than 2 vanish. Finally, when the covariance 
functions (and the connected analogous) are evaluated for 
Xi = X2 = ■■■ = xi v, respectively the moments and the 
cumulants are obtained. 

Continuous random fields and point processes (random 
distributions of discrete points) are both cases of random 
processes. It is possible to relate discrete and continuous 
distributions as we will see in the next paragraphs. 

3.1 Discrete and continuous distributions 

As already mentioned, the theory predicts the statistical 
properties of the continuous matter distribution, while ob¬ 
servations are concerned with the galaxy distribution, which 
is discrete. Different biasing schemes have been introduced 
to relate the (theoretical) distribution of relative mass den¬ 
sity fluctuations to the (observed) distribution of galaxies, 
but the true relation is still unknown. However, indepen¬ 
dently of the presence of bias, it is possible to relate formally 
one to the other. 

In order to recover a continuous distribution from a 
discrete one, define for the discrete process the analogue 
of the density field, the number density field, as a sum over 
Dirac delta functions: 

f(x)=n( x) = h[ 1 + <5(:r)] = y^(5°(x — Xj), 

i 

Y S D {x-Xi)^j, (23) 

where h is the mean number density. 

Conversely, a simple way to generate a discrete distri¬ 
bution from a continuous one is to assume Poisson sampling 
(Peebles 1980, Bertschinger 1992, Fry 1985) by placing, at 
each volume element SV in the sample region, a point with 


probability SP = p(x)5V , where p{x) is the local mean den¬ 
sity. The outcome is a double stochastic process, with one 
level of randomness coming from the random field and the 
second from the Poisson sampling^. It is then possible to 
take into account the effect of discreteness by performing 
a simple substitution in the partition functional or in the 
generating functional for the continuous case. 

In the Poisson model if the probability of finding one 
object in an infinitesimal volume element SV at position x 
is Pi = p(x)SV, then the probability of finding no objects 
is Po = 1 — p(x)5V and the probability of finding more 
than one object is an infinitesimal of higher order. To ob¬ 
tain the probability generating functional for the discrete 
process, we need to generalize this procedure to infinitely 
many small volumes. With this aim, we start from a given 
realization of the stochastic process and divide the volume 
in an infinite set of cells with individual volumes SVi. Then 
we consider the joint probability P{n ( } of finding Ni objects 
in the volume SV i, N 2 in the volume SV 2 , etc., and we define 
the generator of discrete counts in that realization as 

z[J\ = 

{ N t } 

= t 1 - PUhWi + p(*)5%e iJft) ] 
a 

= exp - l)p(xe)5Vz\ . (24) 

£ 

We then let SVt —> 0 in the given realization and average 
over the statistical ensemble of the underlying continuous 
field. This gives the required generating functional for the 
discretized process 



Accounting for the effect of graininess therefore amounts to 
the substitution 


ij -► (e iJ - 1) (26) 

in the generating functional for the continuous process. This 
result is a generalization of the standard procedure leading 
to the generating function for the moments of discrete counts 
in cells (e.g. Fry 1985). The practical utility of this method 
is that it permits calculation of correlated counts in different 
cells, according to the procedures described in Section 3.2 
and Section 3.3 below. 

Actually, there is a subtlety to point out here: the dis¬ 
crete process that we can observe is the one connected with 
the underlying continuous density distribution p(x), but we 
want to study the statistics of the fractional over-density 
field S = p(x)/p — 1 which has several advantages, mainly 
the fact that it has zero mean. To calculate its properties we 
have to modify slightly the above procedure. We first notice 
that the discrete version of the random field 1 + 5 is obtained 
by the substitution 

T This procedure to account for discreteness is a model, and may 
not hold in practice. An example arises in Voronoi tesselations 
(see e.g. Williams, Peacock & Heavens 1991), where objects are 
placed at the intersections of Voronoi planes. This case can even 
result in sub-shot noise power. 
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ij -*■ {e iJ/N - l) N, (27) 

where N is the average number of objects within the sam¬ 
pling volume, in the generating functional of 1 + 5. Next, to 
subtract the terms coming from the mean value 1 , we use the 
fact that the generating functional for the sum of two inde¬ 
pendent random fields is just the product of the individual 
generating functionals. The discrete generating functional 
for 5 is therefore 

Z[J]s = exp i J tf^J’contOz) j Z[J] 1+S (28) 

because the generating functional for a uniform and contin¬ 
uous field of value 1 is simply exp J d 3 xJ CO nt{x)\ • 

In the next section, we give explicit expressions for the 
source function f7, in real and Fourier space, and for contin¬ 
uous and discrete distributions. 


3.2 TV-point correlation functions in real and 
Fourier space 

The quantity that we are dealing with is a complex number 
Z a given by <^<5^(5^ where the vectors form a triangle: 

k\ + /b 2 + k 3 — 0 . 

The average of Z a over the ensemble (of closed trian¬ 
gles) is the bispectrum, but the actual quantity Z a fluctuates 
across the ensemble. In order to quantify this fluctuation it 
would be valuable to estimate its covariance and this in¬ 
volves the quantity (Z a Zp) (see equation |~tGj). 

This quantity is nothing but the six-point correlation 
function (not just the irreducible or connected one) in 
Fourier space (shot noise contribution included) in the par¬ 
ticular case when the six vectors form two triangles: 


VkAk2 .’ ( 29 ) 

where fci + k 2 + k 3 = 0 and k 4 + k 3 + kg = 0 . 

The formalism illustrated in section 3 has been intro¬ 
duced for this reason. Since we want to take into account 
the shot noise contribution, we need to work with the gen¬ 
erating functional for the covariance functions of 1 + S(x) 
which is: 


Z[J] = J V(S)V{5) exp J d 3 xj(x) [1 + 5(x) 
= / exp < i / d 3 xj(x) [1 + 5(x)] 


(30) 


Instead of dealing with functional derivation, we may write 
the ‘external source’ J in a form which allows the func¬ 
tional derivatives to be replaced by normal differentiation 
(Matarrese, Lucchin & Bonometto 1986): 

N 

J(x) SmS D (x ~ Xm). (31) 

m= 1 


and the corresponding connected parts are obtained per¬ 
forming the same kind of differentiation on In Z\J\. 

Smoothing can be incorporated easily. If the smoothed 
field is defined by the convolution: 


^smoothed (x ) 


d 3 x8(x)W (x — x), 


(33) 


the statistics of the smoothed field may be obtained by the 
substitution 


S D (x-Xm) - >W(x-Xm) (34) 

in equation © . For a discrete process the generating func¬ 
tional Z d \J\ is obtained with the substitution (pi]). 

In Fourier space, the Appoint functions may be obtained 
by a straightforward change to the source function, which 
becomes the Fourier transform of the real-space J. If W(k) 
is the transform of the smoothing function, then, for a con¬ 
tinuous field, 

N 

J%(x) = © s m e-*-*W(jL) 

m= 1 

and in the discrete case: 

N 

m=l 

In fact notice that the unsmoothed J'i can be obtained from 
the J d just by substituting, in the exponential, the Dirac 
delta function by its Fourier transform. 

Our ansatz for the generating functional is to assume 
that all correlations vanish above the 3-point function: 


Jt (T) = -iN < exp 



exp 


z\J\ = 

i J d 3 xj(x) - i/ d 3 xd 3 x'J(x)J(x')^ll ln .(x,x') 


— g/ d 3 xd 3 x l d 3 x"j(x)j(x)j(x')^cinn.(x,x,x") 


The validity of this form relies on the fact that we are 
assuming Gaussian initial fluctuations. In the linear regime 
all the connected correlation function of order 3 or higher 
vanish. Allowing then for a quasi-linear evolution, and per¬ 
forming a second-order perturbation theory approximation, 
we allow to become non-zero, but we still rely on all the 
higher-order irreducible correlation functions being negligi¬ 
ble. 

Despite the conceptual simplicity of the algorithm intro¬ 
duced, it is quite cumbersome to perform the 6 th derivative 
of Z d [J], so for practical purposes we use the cluster expan¬ 
sion and express the six-point correlation function in terms 
of the connected parts of lower or equal orders. 


(Si . S 6 ) d/C = 


(5i5 2 ) d clL. (SsSAit n. (5 5 s 6 ) d c lL. 

+ . 

... 15 terms 


+(Si S2) d ^L.{S 3 d i S 5 8 6 ) d c ^. 

+ . 

... 15 terms 


+ VlS2S 3 )± n .{SV5S 6 ) d cinn. 

+ . 

... 10 terms 





(38) 


The covariance functions up to the N th order can then be 
obtained by the following differentiation: 


([1 + 5(*i)] 


[1 + 5(5jv)]) 


d N Z[J] 

dsi....dsN 


Sm= 0 


(32) 


The d/c superscript stands for discrete or continuous case 
and it has been explicitly written to show that it is possible 
to treat the two cases, in configuration or Fourier space, 
on the same footing provided we make the corresponding 
substitution for the external source. 
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For our purposes the calculations have been made di¬ 
rectly in Fourier space for the discrete case giving: 


(SiSj) conn. = (2tt) 3 P(ki) + - S° (hi + kj) , 


(39) 


( SlSmSn)conn. = (27r) 3 8° (kl + km + kn) X 

{ Blmn + ~ [ P{k t ) + P(km) + P(k n )] + J- } . 

In 


(40) 


{SoSpSqSr)conn. = (27r) 3 <5 (k a + kp + k q + k r ) X 

| i [ B( 0+P ) qr + perm. (6 terms)] + 
l. n 

i [( Po+p+q + eye. (4 terms))+ 

Po+P + Po+q + Po+r\ + 

*}■ 


(41) 


(cil . . . <5 6 )conn. = (2n) 3 S D (k 1 + ... + k 6 ) X 

[Si2(3+...+6) + perm. (15 terms) + 

B(i + 2 )(3+4)(5+6) + perm. (15 terms.) + 
Hi(2+3) (4+5+6) + perm. (60 terms)] + 

ttt [Hi + ...-)- Pq + 

n 4 

Pi+2 + perm. (15 terms) + 

Pi+2+3 + perm. (10 terms)] + 



(42) 


The correct indices in the permutations and in the cycli¬ 
cal sequences can also be obtained using graph theory (Fry 
1984, Bertschinger 1992, Matsubara 1995). The expression 
for the six-point discrete correlation function is one of the 
main original results of this work. Its expression is quite com¬ 
plicated, and, it is also important to notice, the discreteness 
terms do not contribute linearly, but they are multiplied 
several times by the signal part (B and P) and other noise 
parts of different order. 


3.3 Selection Function 

The above analysis is valid for volume-limited samples, 
where the mean number density is independent of position. 
In most astronomical catalogues, a selection criterion such 
as a flux limit introduces a position-dependent mean density, 
which adds complication to the analysis. However, it turns 
out that it is very straightforward to include a non-uniform 
selection function in the generating functional approach, by 
simply altering the definition of J. Although we will not 
apply these results in this paper, we present them for com¬ 
pleteness. 

Let the mean density be n(x). We follow Feldman, 
Kaiser & Peacock (1994; hereafter FKP) in defining a fluc¬ 
tuation field 

F(: r) = A w(x)[n(x) — an s (x)], (43) 


where A is a normalisation factor, to be determined later, 
n(x) is the number density of galaxies, n s {x) is the number 
density of a Poisson-sampled synthetic catalogue with mean 
density n(r)/a, and we consider the limit a —> 0, to avoid 
shot noise in the synthetic catalogue. w(x) is an arbitrary 
weighting function, which may be chosen, as in FKP for 
example, to minimise the variance in the estimate of the 
mean power in a shell. 

The Appoint correlations in Fourier space are most eas¬ 
ily calculated by considering the process / = A w(x)n(x) for 
the galaxy field and similarly for the synthetic catalogue 
f s = —a\w(x)n s {x). 

The latter has to be thought as a Poisson sampling of 
a continuous and uniform underlying process whose density 
is diluted by a factor a, therefore the generating functional 
for the process f s is: 


Z a {J) = exp 


and J is given by ([3lJ) or 

The generating functiona 


q 

[exp(— aijfiw) — 1 ]—ax 
a 


(44) 


for the F process is then given 
in terms of the generating functional for the galaxies Z g by 


Z F (J) = Z g (J)Z s (J). 


(45) 


The subtraction of the synthetic catalogue gives us directly 
the Appoint correlations of the zero-mean field F, by the 
methods presented in Section 3.2. Equation ( [l5| ) includes 
(pj|) as a particular case when a —> 0. 


3.3.1 Power spectrum and Bispectrum 


We present here the first two non-trivial terms in the expan¬ 
sion, for illustration, in the case of a volume-limited sample 
and the general case. 

Case (1): Volume-limited catalogue. We choose 



(46) 


so that F(x) = 8(x). In this case the transform of F(x) can 
be used to get an unbiased estimator of the power spectrum: 


(\F K \ 2 ) = V P g (k) + - 


(47) 


where V is the volume of the survey. The 3-point function 
for a closed triangle is 


<*Ci*fc*«3> — V \B g (ki, fe, £ 3 ) 

-t(H 9 1 + Pg 2 + PgS.) + 


(48) 


Comparison with the infinite volume case shows that the 
Dirac delta function is replaced by V/(27r) 3 . Results for dif¬ 
ferent Fourier transform conventions and power spectrum 
and bispectrum definitions are given in the Appendix. 

Case (2): Varying selection function. We choose 

A 2 = I 22 , (49) 


where 



d?xw l (x)ri‘ (x). 


(50) 


Here we find 

(\F%\ 2 ) = PM + ¥-, (51) 
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(-^kl^k2^k3 


> = 


^33 

r 3/2 


B g {kl,ka, fe) + ^-(Pgl + Pg 2 + P 93 ) + ^ 

^33 -/33J 


(52) 


Comparison with the uniform cube shows that to get the 
iV-point functions in the general case the following substi¬ 
tutions may be made: 


Inn _ 1 iN(N-q) 

fNJ2 ’ J NN 


(53) 


4 THE CHOICE OF TRIANGLES IN 

EVALUATING THE BISPECTRUM 

The methods of Section 3 allow us to calculate the means 
and covariance matrix for an arbitrary set of triangles, so 
we may use maximum likelihood to estimate the bias pa¬ 
rameters (equation ^) . We are left with a question of which 
triangles to use as our data. 

The choice of triangles to analyse is of course very wide. 
A likelihood analysis could include a very large number of 
different shapes and sizes, provided that the full covariance 
matrix is calculated. In practice, CPU and memory con¬ 
siderations force one to consider a subset of triangles, and 
the analysis is considerably simplified if the covariance ma¬ 
trix is diagonal. This can be achieved by ensuring that any 
k appears in only one triangle (and, since <5_g = 5 J, one 
must ensure that — k doesn’t appear elsewhere either). This 
throws away some of the information, but the number of 
independent data points is at most 3 times the number one 
uses with this method, so the information content may not 
be reduced significantly. The choice of triangle shape is in¬ 
fluenced by the behaviour of the bispectrum. As we shall 
see, most of the signal comes from the largest k-vectors for 
which second-order perturbation theory holds. We therefore 
make the vectors as large as possible, and for this paper we 
choose only equilateral triangles. As discussed in Section 2, 
a single choice of shape does not allow us to separate 61 and 
62 . Indeed, in the absence of noise, a single triangle shape 
determines the combination &i/(1 + & 2 / 2 J 61 ). In this paper, 
we choose two shapes of triangle, to lift the degeneracy. The 
first is the equilateral configuration, the second we refer to 
as ‘degenerate’, and consists of a repeated vector, and one 
of twice the amplitude and opposite direction. 

In practice the equilateral triangles are obtained as fol¬ 
lows: all the fc-vectors characterized by a polar angle less or 
equal to 60° and a longitude angle less or equal to 180°, are 
considered to be the first vector of a triangle. The second 
and the third are lying symmetrically in a cone of semi¬ 
amplitude 60° centred around —kq. Obviously the choice is 
not unique, and is made at random. The number of actually 
distinct triangles that one can obtain for a given k\ is limited 
by the discreteness of the fc-space and increases with |fc|. A 
similar procedure is followed for the degenerate case, but no 
wavevectors appear in two triangles (of either shape). 

For a discrete grid in a volume-limited cuboidal vol¬ 
ume, neighbouring triangles in fc-space give uncorrelated es¬ 
timates of the bispectrum. This feature was checked by in¬ 
vestigating the values taken by the off-diagonal terms in the 
covariance matrix. 

We have tested the method with an JV-body simulation 
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Figure 2. The power spectrum in the IV-body simulation (solid), 
along with the linear power spectrum (dashed). 


created with Couchman’s (1991) AP 3 M code from a Cold 
Dark Matter (CDM)-like simulation from the Virgo consor¬ 
tium (Jenkins et al. 1996) chosen to match as closely as 
possible the present-day galaxy power spectrum. 

In Fig. 2 we have evidence in the power spectrum for 
the breakdown of linear perturbation theory at A 2 g (k) = 
k 3 P g (k)/(2n 2 ) ~ 1, or at a wavenumber of k ~ 0.3. Fig. 
3 and 4 show the bispectrum (for the galaxy field), in the 
dimensionless form 

x 3 = (^y B kiki ( 54 ) 

and with this convention 

{S 3 ) = [ d cos 9\dtpi ^ k 1 d cos # 2(^2 ^ ^ \ 3 (fci, £ 2 )■ (55) 

J 4-7T 4-7T 

In these figures we see evidence of the breakdown of second- 
order perturbation at k ~ 0.6 for the equilateral triangle, 
and at a much higher k ~ 1 for the degenerate case. The 
non-linear A 2 (fc) ~ 3 and 12 respectively (cf Fry, Melott & 
Shandarin 1995). We use triangles with wavenumbers up to 
k = 0.55 for the equilateral case, and between 0.55 and 1.1 
in the degenerate case in the analysis which follows. 

In practice, we have used the observed power spectrum 
(after shot noise subtraction), rather than the linear power 
spectrum, to calculate the predicted bispectrum (61 = 1 , 
62 = 0), shown in Fig. 3 and 4. This includes second-order 
corrections, but for the leading-order bispectrum term, there 
is no difference. An alternative would be to linearise the 
power spectrum, as in Peacock & Dodds (1994), but this 
requires knowledge of flo- 


4.1 Likelihood analysis 

As a first step we analyze the unbiased distribution of mass 
points, seeking to recover the values 61 = 1 given that 62 = 0 . 

In the case where the covariance matrix is diagonal the 
likelihood function is simplified 
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Figure 3. The measured bispectrum from the simulation (solid), 
with the theoretical prediction from second-order perturbation 
theory shown dashed. Equilateral triangles are used, and no 
wavevector appears in more than one triangle. Error bars are 
errors in the mean for each bin. 


Figure 5. Likelihood of ci = 1/bi where b\ is the linear bias 
parameter, for the (unbiased) ./V-body simulation, for 0.55 < k < 
1 . 1 , using degenerate triangles, and for equilateral triangles (dot¬ 
ted) for k < 0.55. 

by splitting the volume into sub-samples, as is shown in the 
next section. 


£ 

Oxl 



Figure 4. As in Fig. 3, but for ‘degenerate’ triangles k in the 
abscissa is the smaller wavevector. 


£(ci,C 2 ) = 


{2n) M / 2 n a a c 


■ exp < - 


1 \ ' (D a - g 0 

2 Z—d (j2 


(56) 


and very large numbers of triangles can be used. We have 
used as many as possible between the largest wavelengths in 
the cube and the breakdown of second-order perturbation 
theory (determined above). Fig. 5 shows the likelihood for 
fei, assuming that 62 = 0 , from the degenerate and equi¬ 
lateral triangles. The correct value is recovered, within the 
errors, but the error bar in the equilateral case is uncomfort¬ 
ably large. Rather curiously, the error bar can be reduced 


4.2 Expected variance in ci 


We ignore shot noise for simplicity. For equilateral triangles 
we have: 

Ua = (Da) = (Re(5kih2 5 k3)) = c iyP 9 2 <5 D = Cl U ( a- (57) 

The variance in this approximation is 

oi = (Dl) - (Da) 2 = ip 9 3 (<5 D ) 3 . (58) 

The factor of 2 is the same factor as in the transition from 
complex Z a to real D a in equation (^|). Note that of the 
fifteen signal terms in the covariance matrix (Z a Zp), only 
one survives. For a discrete transform in a cube of volume 
V, the Dirac delta function is replaced by V/(2n) 3 , as in 

©■ 

Since cr 2 is independent of ci, the error on ci is 


d 2 In £ 


dc\ 


= £ 


\n 


(l)l2 


(59) 


In practice we are dealing with more than a single mode, 
and moreover, in fc-space the modes form a discrete set, with 
density of states g = V/(2 / k) 3 . The number of uncorrelated 
triangles in a thin shell of width 8(Ink) is: 

^4nk 3 5(lnk)g^, (60) 

where the factor i arises from the reality of the density 
fluctuation field and the factor | from the requirement of 
uncorrelated data (no wavevectors in common between any 
two triangles). 

Considering contributions from all the shells, and mov¬ 
ing to the continuous limit this becomes: 
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Figure 6. Likelihood of c\ = l/b\ for several sub-volumes, each 
containing 1/8 the volume, for degenerate triangles for the same 
AT-body simulation as Fig. 4. 


Figure 7. Joint likelihood of c\ = l/b± and C 2 = 62 /^ 1 * Con¬ 
tours contain 68.3 and 90 percent of the a posteriori probabil¬ 
ity assuming the likelihood is a bivariate gaussian, and the prior 
probabilities for c\ and C 2 are uniform. 


-2 


cr 


ci 


V 

12t r 2 



d( In k) 


[Mb > (1) ] 2 ,_ 3 

cr 2 (fc) ’ 


(61) 


where fc ma x is set by the breakdown of perturbation theory, 
and fc m in is set by the size of the sample. 

Inserting the appropriate expressions in the previous 
equation one obtains: 


ffei2 ' = ‘'§ ff “ (fcm “)’ ( 62 ) 
where 


&0 (^max) — 2^2 


d(ln k)P g k 3 


(63) 


is the variance in the galaxy field for a top hat filter in 
fc-space truncated at fc ma x- It is also truncated at the sam¬ 
ple size at low-fc, which makes the definition slightly non¬ 
standard. Hence 


a 


ci 


1.01 

n’o(^-max) 


(64) 


which gives the error of 1.48 seen in Fig. 5, for the cutoff at 
k = 0.55. A curious feature of this analysis is that the error 
is independent of the size of the survey, motivating the split¬ 
ting of the volume into independent smaller units. This loses 
some long-wavelengtli modes, but these contain little infor¬ 
mation anyway. Fig. 6 shows the likelihood curves for several 
sub-volumes, and comparison with Fig. 5 demonstrates that 
the error is independent of the size of the subvolume. This 
is somewhat reminiscent of the equally counter-intuitive re¬ 
sult that the variance in an estimate of the power spectrum 
is also independent of the size of the survey, motivating a 
splitting of the volume as a possible way to increase signal- 
to-noise in that case (Press et al. 1992). 


4.3 Determining bi and ?>2 

For simplicity, the error analysis of the previous section as¬ 
sumed that the quadratic bias term 62 is zero. In practice, 
this may not be true, and we need to estimate both 61 and 
62 from the data. To lift the degeneracy between 61 and 62 , 
we need to consider at least two triangle configurations. The 
optimal configurations are not easy to determine, and will 
certainly depend on the power spectrum, as anything other 
than the equilateral configuration involves wavenumbers of 
different magnitude. In addition, each triangle shape leads 
to approximate determination of a degenerate combination 
of bi and & 2 , so ideally we would like these degenerate lines 
to cross at large angle. Adding the degenerate triangles to 
the equilateral configuration is particularly effective, as each 
degenerate triangle removes only two wavevectors from the 
available set, and also second-order perturbation theory ap¬ 
pears to work well into the nonlinear regime. Fig. 7 shows 
the joint likelihood for ci and C 2 for the limits on k pre¬ 
viously stated. The true solution is nicely within the error 
bars. 


5 DISCUSSION 

We have presented here the first steps towards the goal of 
using the bispectrum to measure the linear bias parameter b 
from large-scale structure data. The aim of this is ultimately 
to combine measurements of b with estimates of f3 = Q.q' 6 /b 
from redshift distortion studies to get an unambiguous es¬ 
timate of the density parameter flo- In this paper we have 
presented a method for calculating the covariance matrix for 
bispectrum estimates, thus allowing maximum likelihood to 
be used to estimate b, and tested the methods on N-body 
simulations. The noise calculation is a general method which 
will be useful for calculating the properties of any high-order 
statistics, and is the major new advance of this paper. We 
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have also shown how the signal-to-noise of the bias estimate 
may be improved by splitting the volume into subsamples. 

It is worth making some remarks about the shape infor¬ 
mation, as this is crucial in lifting a degeneracy which exists 
physically between the effect of gravitational evolution and 
nonlinear bias (in the analysis, this manifests itself as a de¬ 
generacy between the linear and quadratic bias terms). In 
essence, we are exploiting the fact that gravitational insta¬ 
bility skews the density field away from its initial Gaussian 
symmetry. However, this can also be achieved by a nonlinear 
bias, so we need a further property to distinguish the two. 
It is the shape information which allows this. Although we 
work in Fourier space, the argument is most easily stated in 
real space. The effect of a nonlinear bias, operating on the 
Gaussian field, is to shift iso-density contours up (or down), 
but maintaining the shape of the contour. Gravitational evo¬ 
lution, however, will change the shapes, usually leading to 
flattening of collapsing structures (e.g. Lin, Mestel & Shu 
1965). By utilising shape information, we can, in principle, 
decide to what extent nonlinear biasing and gravitational 
evolution are responsible for the skewed density field. 

A number of unresolved issues remain, and these will 
be addressed in a subsequent paper. There are two ma¬ 
jor ones: the practical effects of a varying selection func¬ 
tion, and the effect of redshift-space distortions, arising from 
the fact that the redshift is an imperfect distance indicator 
because of the contaminating effects of peculiar velocities. 
This effect can be regarded as being split into two compo¬ 
nents: a large-scale distortion resulting from coherent inflow 
into overdense regions, and the ‘Fingers-of-God’ arising from 
virialised, highly non-linear structures. Since we have argued 
that it is effective to use subsamples, we would expect to be 
able to use Kaiser’s (1987) distant-observer approximation, 
for a deep survey, to model the large-scale distortions (see 
also Hivon et al. 1995). Of more concern is that most of the 
signal used to measure the bias parameters comes from the 
mildly non-linear regime, and it may be that simple models 
for the Fingers-of-God (such as the incoherent velocity dis¬ 
persion introduced by Peacock (1992)) may not be adequate 
for our purposes. 

Assuming these problems can be dealt with, the accu¬ 
racy of the determination of b will be dependent on the num¬ 
ber of subvolumes we make, and this is determined by how 
small the subvolumes may be. From the numerical simula¬ 
tions, we find that the subvolumes degrade significantly the 
bispectrum estimate if the size of the box is less than about 
the nonlinear wavelength 27 t/A:nl, or about 20 h, -1 Mpc for 
IRAS galaxy data. One would therefore expect rather weak 
constraints on b from the IR AS PSCz survey, but a much im¬ 
proved error of under 10% from the Sloan digital sky survey 
or Anglo-Australian 2 degree-field survey. These are rough 
estimates, based on the shot noise and depth of each survey; 
it will require more detailed study to see if this accuracy can 
be achieved or surpassed. 
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APPENDIX A: FOURIER TRANSFORM 
CONVENTIONS 


There are many possibilities for placing of factors of 2-7T etc. 
in the definitions of the power spectrum, bispectrum and 
Fourier transform. These are characterised by the constants 
A, D and E below. Our choice is A — 1, D E = (2-7t) 3 , 
which seems to yield the minimum complexity in the for¬ 
mulae. With the addition of E, this is the same conven¬ 
tion as FKP, but differs from Bertschinger (1992), who uses 
A = (2 n)~ 3 , D = E= 1 . 

= A J d 3 xS(x) exp (—ik.x) 

<<^A 2 > = DP(k)S D (ki + k 2 ) 

(kikjts) = EB(k 1 ,k 2 ,k 3 ) 5 D (k 1 + k 2 + k 3 ). 

(Al) 


We present below some of the simpler results, including shot 
noise, with general Fourier conventions for workers whose 
convention is different from ours. For the case of a varying 
selection function, the power is 


(\ F k\ 2 ) = > 2 


P(k)DI 22 

(2tt) 3 


+ A I 2 \ 


and the 3-point function is 


(A2) 


( F kl F k2 F ks) 


A 3 


EI33 

J2FY 


B(ki, k 2 , £3)+ 


AD/32 

(2tt) 3 


(Pi + P2+ Pi) + A 3 131 


(A3) 


The conventions may be related to the inverses 

P(k ) = J ( j3 x ^(2) exp(—ik .x) 

B(k!,k 2 ,k 3 ) = (27r ^ 3 J d 3 x 12 d 3 x 23 

V 3) (■ x , x + xi 2 ,x + * 23 ) exp(— ik 2 .xi 2 — ik3-X\3) 

(A4) 

(the latter for a closed triangle). The source function for 
the generating functional in Fourier space is modified in the 
continuous case from (|35| ) to: 

N 

J%(x) = A ^ Sm.e~^' x W (km)- (A5) 
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